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CHANNEL  MAINTENANCE  BY  TRAINING  STRUCTURES:  GUIDANCE 


FOR  NUMERICAL  MODEL  MESH  DEVELOPMENT 


PART  I:  INTRODUCTION 

Background 

1.  Natural  rivers  and  estuaries  are  quite  shallow  relative  to  their 
width;  typical  depth-to-width  ratios  are  less  than  1/10,  which  implies  that 
the  flow  distribution  is  dominated  by  bottom  friction  rather  than  the  drag 
from  the  channel  banks  (Chow  1959).  Man-made  training  structures,  however, 
are  a  sufficient  anomaly  in  this  environment  to  reverse  the  relative  impor¬ 
tance  of  these  two.  The  training  structure  forces  large  velocity  and  water- 
surface  gradients  and,  perhaps  more  importantly,  severe  changes  in  these  gra¬ 
dients.  Thus  the  structure  will  profoundly  influence  the  flow  pattern  over  a 
large  region. 

2.  The  complexity  of  the  flow  fields  produced  by  training  structures 
often  precludes  analyses  short  of  physical  or  numerical  modeling.  Two- 
dimensional  numerical  models  founded  upon  the  shallow-water  equations  are 
common  and  are  perhaps  the  simplest  model  capable  of  reasonably  examining 
training  structures.  This  report  discusses  the  limitations  of  the  numerical 
representation  of  these  equations  on  applications  with  training  structures. 

3.  The  focus  is  upon  the  mesh  development.  Too  often  the  factors  in 
forming  an  adequate  mesh  are  considered  only  after  precision  or  stability 
problems  are  encountered,  perhaps  in  final  runs.  A  properly  fabricated  mesh 
from  the  onset  will  save  time  and  expense  avoiding  costly  reruns  in  a  piece¬ 
meal  fashion. 


Objectives 

4.  The  objectives  of  this  investigation  are  to  describe  the  influence 
of  geometry  and  the  mesh  upon  the  precision  of  shallow-water  numerical  model 
results  and  to  produce  simple  instructions  and  rules  of  thumb  for  mesh  devel¬ 
opment.  With  this  in  mind,  the  following  elements  of  mesh  production  are 
considered; 
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a.  Resolution. 

b.  Mesh  orientation. 

c.  Skewness  in  the  grid. 

d.  Oscillation  suppression. 

e.  Bathymetric  effects. 


iiroach 


5.  Assessments  are  made  from  a  set  of  first-order  linear  model  equa¬ 
tions  applied  in  a  Galerkin  finite  element  framework.  All  calculations  are 
for  steady  model  equations;  therefore,  no  discussion  of  time  discretization 
errors  such  as  numerical  dispersion  is  included.  Steady  models  provide  an 
uncomplicated  method  to  evaluate  grid-related  problems.  The  first  three  items 
listed  in  paragraph  4  are  resolved  via  the  calculation  of  the  truncation  error 
in  a  one-  or  two-dimensional  finite  element  model  using  a  linear  basis.  The 
last  two  items  are  determined  from  the  amplification  caused  by  a  one¬ 
dimensional  finite  element  representation  of  a  steady-state  linearized 
shallow-water  model. 
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PART  II:  NUMERICAL  BACKGROUND 


The  Shallow-Water  Equations 


6.  The  shallow-water  equations  in  one  dimension  may  be  written 


H,.  +  UHjj  +  HUX  =  0 


(1) 


Ut  +  UUX  =  i/U^  -  g(H  ♦  x)x  (2) 

where* 

H  =  depth 

t  =  time 

U  =  velocity 

x  =  distance 

v  =  apparent  viscosity 

g  =  acceleration  due  to  gravity 

z  =  channel  bottom  elevation 

and  the  subscripts  indicate  differentiation.  The  basic  assumptions  in  this 
formulation  are  mild  bed  slope  and  hydrostatic  pressure  distribution.  The 
subsequent  analysis  is  for  steady  flow  conditions,  which  can  be  obtained  from 
the  previous  equations  by  simply  dropping  the  terms  Ht  and  Ut  .  These 
equations  can  be  derived  from  Euler's  equations  in  a  systematic  fashion  as  in 
Stoker  (1957)  . 


Numerical  Considerations 

7.  The  numerical  representation  of  the  shallow-water  equations  can  be 
subject  to  significant  discretization  errors,  the  introduction  of  which  is  a 
concern.  The  analytic  solution  of  the  differential  equations  (shallow-water 
equations)  preserves  the  relationship  among  derivatives  of  variables  specified 
by  differential  equations  at  EVERY  point  in  the  domain.  The  numerical  model, 
on  the  other  hand,  is  solved  based  upon  a  finite  number  of  points.  If  the 


*  For  convenience,  symbols  and  unusual  abbreviations  are  listed  and  defined 
in  the  Notation  (Appendix  B) . 


model  is  properly  formulated,  then  as  the  number  of  mesh  points  increases  to 
infinity,  the  analytic  solution  is  recovered.  This  property  is  termed 
consistency . 

8.  As  a  result  of  the  limited  amount  of  information  involved  in  a 
single  equation,  truncation  error  is  introduced.  In  the  case  of  a  global 
expansion  method,  all  data  in  the  domain  are  used  for  each  numerical  equation, 
thus  providing  a  low  truncation  error.  In  finite  difference  and  finite  ele¬ 
ment  methods,  only  a  limited  number  of  neighboring  points  are  used,  reducing 
computational  expense  but  producing  much  more  truncation  error.  The  trunca¬ 
tion  error  associated  with  a  particular  numerical  scheme  is  found  by  a  Taylor 
series  expansion  of  each  nodal  value  about  a  common  node  location.  In  the 
following  work,  simple  model  equations  are  used  to  assess  the  truncation  error 
in  two  spatial  dimensions  for  linear  basis  using  triangular  and  quadrilateral 
finite  elements.  The  extension  to  finite  difference  is  readily  apparent.  The 
following  sections  shall  consider  the  truncation  error  of  various  finite  ele¬ 
ment  patches  for  a  simple  first-order  model  equation.  Stability  problems  that 
may  arise  as  a  result  of  resolution  and  by  the  variation  of  bathymetry  will 
then  be  discussed. 


Background  for  Finite  Elements 

9.  The  Galerkin  finite  element  technique  as  applied  to  differential 
equations  will  be  approached  through  the  idea  of  the  approximation  or  interpo¬ 
lation  of  a  function  by  some  simpler  set  of  functions.  This  in  turn  leads 
directly  to  application  to  differential  equations  in  which  relationships 
between  the  function  and  its  derivatives  are  to  be  maintained.  a  triable 
such  as  concentration  or  velocity  components  can  be  approximated  as  a  combina¬ 
tion  of  functions.  For  example,  the  concentration  C(x)  can  be  approximated 
as  a  linear  function  as 

C(x)  =  a:x  +  ,  for  0  <  x  <  1  (3) 

where  ac  and  aj  are  constants.  This  can  be  rewritten  equivalently  as 

C(x)  =  Cc(l  -  x)  +  C:x  .for  0  <  x  <  1  <*) 
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where 

C0  =  C(o) 

C  i  =  C(l) 

are  values  at  x  =  0  and  x  =  1  ,  nodal  values.  Equation  4  can  then  be 
simplified  to  the  form 


C(x) 


2 


=  E 


(5) 


where 

<Mx)  -  (1  -  x) 

02(x)  =  x 

j  =  an  index 

It  should  he  apparent  that  if  a  higher  order  approximation  is  desired,  it  is 
necessary  only  to  include  additional  functions  and  corresponding  nodal  values. 
The  particular  type  of  function  in  Equation  5  is  called  a  Lagrange  polynomial; 
other  function  types  can  be  developed,  but  these  equations  are  in  common  use 
and  will  be  used  exclusively  throughout  this  report.  The  functions  used  to 
approximate  the  variable  C(x)  are  termed  interpolation  or,  alternatively, 
basis  functions. 

10.  The  finite  element  approach  is  to  apply  these  interpolation  func¬ 
tions  over  each  element  or  region  that  subdivides  the  total  region.  In  this 
manner  the  global  solution  can  be  approximated  precisely  without  the  need  for 
a  high-order  polynomial.  The  total  basis  function  for  a  node  will  then  be  the 
combination  of  functions  for  this  node  in  all  elements  joining  this  node. 

11.  As  an  example  of  an  application,  consider  the  approximation  of  the 
function  sin  (x)  over  the  interval  0  <  x  <  *  by  linear  polynomials.  Only 
two  elements  will  be  used,  though  any  number  could  be  used.  The  approximation 
may  be  written  as 


C(x) 


(6) 


where 
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X  “  Xr 

C  ( x )  =  _ — -  ,  for  Xl  <  x  < 

Xr  XL 


(7) 


Within  each  element,  L  and  R  denote  the  left  and  right  node  locations. 
The  basis  functions  within  the  element  are: 


MO  =  l  -  e 
MO  =  I 


where  <f>  is  the  one-dimensional  linear  shape  function.  To  form  the  complete 
basis  function  <p1  ,  for  example,  combine  <f>R  from  element  1  with  4>l  from 
element  2  (Figure  1) . 

Node  0  l  ? 


Figure  1.  Finite  element  grid  for  approximation  of  sin  x 


12.  Further  assume  that  CQ  =  C2  =  0  ,  which  precisely  matches  the 
function  sin  (x)  at  the  boundaries  and  simplifies  the  problem  of  finding  the 
proper  value  of  Cx  .  An  obvious  way  to  choose  C]^  is  again  to  simply  let 
Cj  =  sin  ( 7t/2 )  =  1.  The  resulting  approximation  would  be  as  shown  in 
Figure  2. 


-1  12  3  4 


x 


Figure  2.  Direct  interpolation 
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As  one  can  see,  chis  representation  leaves  something  to  be  desired.  It  is 
consistently  low.  Another  method  would  be  to  define  error  as  e  =  [f(x) 

-  ^^(xlCj],  where  f(x)  is  the  function  to  be  approximated  (sin  (x)  in  this 
case)  and  £<^(x)C.j  is  the  approximate  function.  Minimize  the  error  norm 

3 

given  by 

li e  1  =  Jn  [f(x>  -  £  <f>i  (x)Cj  ]2  dx  (8) 

where  fl  is  the  domain,  by  taking  the  partial  derivative  with  respect  to  each 
unknown  Cj  .  The  result  is 

f  4>  i(x)  [f(x)  -  £  <^>3  (xlC^  ]  dx  =  0  ,  for  each  i  (9) 

Jn  j 

which  is  a  set  of  algebraic  equations  to  solve.  An  alternative  way  to  write 
this  is 

,  e)  =  0  (10) 

where  this  is  termed  the  inner  product  form.  This  indicates  that  the  error  is 
made  orthogonal  to  the  functions  4>i  termed  the  test  functions.  Here  it  is 
shown  that  if  the  test  functions  are  chosen  to  be  the  same  as  the  interpola¬ 
tion  or  basis  functions,  the  error  is  minimized.  The  result  for  the  example 
is  Cx  -  1.216.  Figure  3  illustrates  that  this  is  a  better  choice  than  the 
previous  case. 


X 

Figure  3.  Galerkin  approximation 
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The  concept  of  using  the  same  test  functions  as  basis  functions  is  the 
Galerkin  technique  in  finite  elements. 

13.  This  technique  needs  to  be  applicable  to  differential  equations. 
For  example,  consider  the  one-dimensional  (1-D)  steady-state  convection 
diffusion  equation 


uc*  -  DCxx  =  0  ,  for  0  <  x  <  1  (11) 

where  D  is  the  diffusion  coefficient.  Note  that  here  the  differential  equa¬ 
tion  explicitly  states  a  relationship  between  the  function  and  its  various 
derivatives.  In  this  case  there  is  a  relationship  between  the  function  slope 
and  curvature.  The  approximation  needs  to  preserve  this  relationship  as 
nearly  as  possible.  The  straightforward  application  is  to  let  C(x) 

=  ^(/^(x)^  .  The  Galerkin  finite  element  approximation  is  then  written  as 

3 


[*i  .  U  £  ^(X)C.  -  D  £  //Cj]  =  0  ,  for  each  <f>.  (12) 

j  3 

where  the  prime  indicates  differentiation.  TMs  notation  avoids  subscripts  j 
and  x  for  different  operations.  It  would  appear  that  <£(x)  must  have  at 
Least  second  derivatives  that  are  nonzero.  However,  <£(x)  can  be  linear  by 
using  integration  by  parts  putting  a  derivative  on  the  test  function.  It 
should  be  noted  that  it  can  be  shown  that  the  Galerkin  technique  provides  an 
"optimum"  approximation  when  dealing  with  the  function  or  any  even-numbered 
derivatives.  Unfortunately,  in  the  case  of  the  convection  operator,  which 
involves  first  derivatives  only,  there  is  no  "norm"  to  be  minimized  by  using 
the  Galerkin  technique.  However,  this  technique  is  commonly  used  for  convec¬ 
tion  problems  as  is  the  case  in  TABS-2 . 

14.  The  terms  have  been  supplied  that  are  necessary  to  understand  the 
finite  element  examples,  which  follow. 
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PART  III:  RESULTS  AND  DISCUSSION 


15.  The  steady— state  shallow-water  equations  are  often  used  to  repre¬ 
sent  estuarine  or  riverine  conditions  in  which  convection  is  much  more  impor¬ 
tant  than  diffusion.  Therefore  a  simple  model  equation  that  is  only  first 
order  shall  be  considered: 

4-o  <13> 

This  function  can  be  thought  of  as  convection  of  some  species  concentration 
f  ,  in  which  the  convection  velocity  has  a  value  of  1  in  the  x-direction  only. 
This  example  model  equation  will  be  used  to  evaluate  1-D  and  two-dimensional 
(2-D)  model  schemes. 


Resolution 

16.  The  resolution  supported  in  a  region  determines  the  precision  of 
the  numerical  results.  However,  there  are  few  definitive  suggestions  about 
the  grid  size  and  the  rate  of  expansion.  In  addressing  the  expansion  rate 
question,  consider  the  truncation  error  for  the  expansion  of  a  1-D  finite 
element  representation  of  fx  =  0  .  Representing  this  numerical  equation  in 
inner  product  form  gives 


(*i(x)  .  £^j(x)f.,)  =  0  ’  for  a11  i 


(14) 


The  notation  simply  means  the  integral  over  the  domain  of  the  product  of  the 
two  functions  (separated  by  the  comma).  The  function  f(x)  is  represented  by 
an  expansion  in  the  basis  function  ^(x)  ,  i.e.,  f(x)  -  £  ^(x)fj  ,  and  the 
prime  indicates  differentiation.  The  subscript  indicates  the  nodal  index,  as 
shown  in  Figure  4.  The  inner  product  uses  the  same  test  function  as  that  used 
as  a  basis  for  the  function  f(x)  ,  thus  a  Galerkin  weighted-residual 
statement . 

Li  L2  L3 

•  •  "•  ■■■  '  ♦— - 

0  1  2  3 

1-D  GRID 

Figure  4.  1-D  grid 
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17.  If  it  is  considered  that  the  gradients  in  water  levels  and  veloci¬ 
ties  and  changes  in  gradients  might  be  quite  large  at  a  dike,  then  the  small¬ 
est  elements  need  to  be  there.  However,  not  being  able  to  afford  this  resolu¬ 
tion  throughout  the  domain,  the  mesh  size  needs  to  be  expanded  in  a  prudent 
fashion  such  that  the  error  would  not  be  amplified  severely.  With  this  in 
mind,  the  truncation  error  associated  with  the  boundary  and  at  a  typical  loca¬ 
tion  within  the  field  shall  now  be  considered.  A  single  algebraic  equation 
results  from  each  selection  of  the  test  function  ^A(x)  .  The  function  is 
nonzero  only  over  the  elements  that  contain  node  i  .  Thus,  the  equation 
includes  only  a  few  nodal  values. 

18.  If  i  -  0  is  chosen,  then  the  algebraic  equation  associated  with 
the  boundary  is  selected.  As  implied  by  Figure  4,  the  basis  chosen  is  linear 
with  the  resulting  algebraic  equation  of 


f,  -  f0  =  0  0-5) 

Expanding  this  expression  in  a  Taylor  series  expanded  about  x0  and  solving 
for  fx  gives  the  truncation  error 


Truncation  error  =  -  1/2  Lj  -  1/6  hf  f„„  -  .  .  .  (16) 


where  L  is  the  distance  between  adjacent  nodes.  Performing  the  same  opera¬ 
tions  on  a  typical  test  function  ^  gives  the  algebraic  equation  and  the 
truncation  error 

fi-i  -  fi-i  =  0  (17) 


Truncation  error 


-1/2(1^  -  1*)^  -  1/6 


Lj.i  »  Lj3 

1*1  *i  +  1* 


fxxx 


(18) 


Therefore,  assuming  f**  at  the  boundary  is  greater  than  or  equal  to  that  out 
in  the  interior,  if  the  element  size  varies  as  -  iLx  (the  element  size 

increases  such  that  the  second  element  is  length  21^  ,  the  third  3L:  ,  and 
so  on),  the  precision  in  the  mesh  is  no  worse  than  that  of  the  boundary. 
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Orientation  and  Skewness 


19.  Again  considering  the  truncation  error  associated  with  a  simple 
model  problem  fx  =  0  ,  the  2-D  finite  element  representation  is  investigated 
as  either  a  set  of  triangles  or  quadrilaterals  in  the  vicinity  of  a  hypotheti¬ 
cal  dike.  The  simple  model  equation  is  analogous  to  the  steady-state  trans¬ 
port  of  f  when  the  velocity  is  along  the  x-direction.  This  is  completely 
general  in  that  the  coordinate  system  could  be  rotated  to  match  the  flow 
direction  and  result  in  this  equation.  Figure  5  illustrates  this  case. 


The  flow  is  in  the  x-direction  and  the  dike  is  perpendicular  to  the  flow, 
i.e.,  in  the  y-direction.  Note  that  downstream  of  the  dike  tip  very  large 
derivatives  are  expected  in  the  y-direction  but  not  in  the  x-direction.  Thus 
a  truncation  error  in  which  only  y-derivatives  are  present  will  be  emphasized. 

20.  The  finite  element  statement  and  the  region  or  patches  over  which 
it  will  be  applied  are  as  follows: 

(^(x.y)  ,  £  ^  (x,y)  f^)  =  0  ,  for  each  i  (19) 


where 

4>  -  2-D  linear  basis  functions,  either  triangle  or  quadrilateral 
elements 

'  -  partial  derivative  with  respect  to  x 
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Figure  6.  Typical  2— D  finite  element  patches 


The  patches  shown  in  Figure  6  result  from  the  nonzero  region  of  the  test  func¬ 
tion  <f> 5  .  This  is  the  only  region  that  needs  to  be  considered  since  in  the 
finite  element  method  when  i  —  5  ,  the  product  of  this  test  function  is  being 
integrated  with  the  partial  differential  equation  and  so  the  only  contribution 
to  the  resulting  algebraic  equation  is  from  the  nonzero  region  of  the  test 
function  <j>$  ,  as  shown  in  the  following  equations: 
a.  Triangles : 


(f6  -  f*)(Ya  -  Y2)  +  <f2  -  fB)(Y6  -  YO  -  0  (20) 

b .  Quadrilaterals : 

fi(Y2  -  Y*)  +  f2(-Yx  +  Y3  -  Y4  +  Yg)  +  f3(-Y2  +  Y6)  +  f*(Yj 
+  Y2  -  Y7  -  Y8)  +  f6(-Y2  -  Y3  +  Ye  +  Y9)  +  f7(Y4  -  Ye) 

+  f,(Y*  -  Yg  +  Y7  -  Y9)  +  f9(— Yg  +  Y8)  =  0  (21) 


where  Y  is  distance  normal  to  the  x-axis. 

21.  These  algebraic  equations  can  be  used  to  calculate  the  truncation 
error  by  expanding  in  a  Taylor  series  about  the  location  of  node  5.  It  is 
assumed  for  simplicity  that  the  metrics  involved  in  the  geometric  transforma¬ 
tion  are  exact.  In  this  manner,  Figure  7  was  developed  showing  the  truncation 
error  for  the  triangular-  and  quadrilateral-based  finite  element  patches  under 
three  circumstances. 

22.  Figure  7a  indicates  the  truncation  error  for  the  case  in  which  both 
grids  form  an  orthogonal  patch  oriented  with  the  flow.  The  truncation  error 
of  the  quadrilateral  patch  includes  an  additional  term  beyond  that  of  the 
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•1/6  L2^ 


•1/6  L2(f)o(x  +  *xyy^ 


a.  Oriented  with  the  flow 


■1/6  L  (fyyy  +  Sccf^^y  •  Gtfyyy) 


-1/6  L2(1mx  ♦ 1^) 


b.  Possible  rotation 


0.5a 

1 

0.5a 

1 

1 

<£>  — - 

1 

EE? — 

1/6  L  ^XXX  +  ^xxy  ■  afyyy) 

-1/6  *-2(fxxx  +  2afXXy  +  fjfyy) 

c.  Skewness 

Figure  7.  Leading  terms  of  truncation  error  for  fx  -  0 
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triangular  mesh.  It  is  not  a  term  fyyy  ;  thus  the  truncation  error  is 
probably  small  for  both  cases  (remember  that  gradients  in  the  x-direction  are 
small).  Possible  rotation  is  shown  in  Figure  7b,  and  skewness  (collapse  of 
the  element)  in  the  patch  is  shown  in  Figure  7c.  Both  cases  deal  with  a  small 
rotation  or  skewness  angle  of  a  ,  where  0  <  a  «  1  . 

23.  The  inclusion  of  a  slight  rotation  results  in  two  additional  terms 
in  the  truncation  error  of  the  triangular  case.  More  importantly,  one  of 
these  terms  involves  fyyy  ,  which  could  be  quite  large  in  the  vicinity  of  the 
dike  tip.  The  quadrilateral  patch,  on  the  other  hand,  had  no  change  in  trun¬ 
cation  error  for  this  condition. 

24.  The  effect  of  a  small  degree  of  skewness  was  identical  to  the  small 
rotation  for  the  triangular  patch,  including  the  potentially  large  term 

fyyy  .  The  quadrilateral  patch  with  this  slight  skewness  was  changed  to 
include  an  additional  term  fMy  . 

25.  It  may  then  be  concluded  that  for  best  results,  gi.ids  in  the  vicin¬ 
ity  of  the  dike  tip  should  be  aligned  with  the  flow  direction  and  orthogonal. 
If  this  is  not  possible,  it  would  be  preferable  to  use  quadrilateral  elements 
in  this  immediate  vicinity. 


Oscillation  Suppression 

26.  A  common  problem  with  shallow-water  numerical  models  is  the  appear¬ 
ance  of  spatial  oscillations  around  steep  changes  in  the  flow  variables.  This 
is  actually  a  problem  in  transport  models  as  well  and  results  from  the  first- 
order  derivatives  in  the  differential  equations.  To  illustrate  the  problem, 
consider  the  case  of  the  steady  convection-diffusion  equation 


Uf,  -  Dfjoj 


(22) 


where 

U  -  convection  velocity 
x  -  coordinate  direction 

Subscripts  indicate  differentiation  with  respect  to  x  .  Representing  f  as 
discrete  nodal  values  in  a  linear  basis  gives 
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(23) 


f(x)  = 

j 

where  4>  is  the  interpolation  function.  The  algebraic  equation  set  resulting 
from  the  Galerkin  finite  element  representation  of  the  convection  diffusion 
equation  for  a  linear  basis  and  a  uniform  grid  is 


<fi*l  -  f3-l)  +  |  <*1.1  -  2fi  +  *!-!>  =  0  <24> 

where  R  —  UL/D  and  is  known  as  the  Peclet  number  or  the  cell  Reynolds 
number.  If  it  is  assumed  that 


fj  =  c  pi  (25) 

where 

c  =  a  constant 
p  =  the  numerical  root 

there  are  two  roots  p  =  1  ,  (2/R  -  l)/(2/R  +  1)  .  The  first  root  implies 
that  f  is  a  constant;  the  second,  however,  allows  variability  in  f  .  But 
note  that  if  p  is  negative,  the  model  representation  will  oscillate.  Thus 
it  must  be  noted  what  restrictions  guarantee  that  p  will  be  positive. 
Limiting  the  minimum  value  to  zero,  the  relationship  to  prevent  oscillation  in 
the  model  is 


R  -  UL/D  <  2  (26) 

Therefore,  it  is  possible  to  refine  the  grid  to  prevent  oscillations,  provided 
the  diffusion  is  not  negligible. 

27.  The  extension  to  the  shallow-water  equations  is  not  straightforward 
as  there  is  a  set  of  equations,  one  of  which  (conservation  of  mass)  contains 
no  diffusion  term  at  all.  Therefore,  there  is  no  means  to  damp  water-surface 
oscillations  in  a  model  using  linear  basis  for  both  water-surface  and  veloc¬ 
ity.  It  has  been  shown  (Platzman  1978,  Walters  and  Carey  1983)  that  by  using 
a  mixed  interpolation  (quadratic  for  velocity  and  linear  for  water  surface) 
the  2L  wavelength  oscillation  does  not  appear  in  the  water  surface.  This 
representation  is  used  in  many  shallow-water  models,  in  particular,  the  TABS-2 
hydrodynamic  model  RMA-2V  (Thomas  and  McAnally  1985).  The  velocity 
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oscillation  is  controllable  via  the  diffusion  and  mesh  resolution.  The  rela¬ 
tionship  for  small  Froude  number,  F  =  U/(gH)1/z  ,  reverts  to  the  same  rela¬ 
tionship  found  for  the  convection-diffusion  model  (for  details  see  Appendix  A) 


R  =  _  <  2  (27) 

u 

where  L  is  the  distance  between  velocity  nodes.  In  the  2-D  case,  it  is 
important  to  note  that  the  streamwise  diffusion  or  viscosity  and  resolution 
are  the  important  parameters  in  suppression  of  oscillations.  It  should  also 
be  noted  that  while  the  2L  wavelength  oscillation  in  the  water  surface  will 
be  suppressed,  longer  wavelength  oscillations  can  occur. 

Bathymetric  Effects 

28.  In  model  applications,  numerical  model  behavior  appears  to  be  much 
less  stable  when  the  bathymetry  changes  rapidly.  Limitations  are  imposed  by 
the  model  equations  themselves  due  to  the  mild  slope  assumption  and  the  lack 
of  variation  of  velocity  over  depth.  However,  this  instability  is  in  fact 
numerically  induced. 

29.  In  addition  to  the  oscillation  suppression  by  diffusion  and  resolu¬ 
tion  discussed  previously,  which  were  evaluated  for  flat  bottom  conditions,  a 
variable  bathymetry  must  now  be  considered.  The  problem  will  be  simplified  by 
excluding  viscosity  and  using  linear  basis  for  both  water  surface  and  veloc¬ 
ity,  the  reasoning  being  that  the  purpose  of  this  exercise  is  illustrative, 
and  can  be  accomplished  even  without  these  additional  complications. 

30.  The  1-D  steady-state  shallow-water  equations  may  be  written  as: 

UHX  +  HUX  -  0  (28) 

gHx  +  UUX  -  -gzx  (29) 

where  U  is  the  longitudinal  velocity.  These  equations  can  be  linearized  by 
considering  small  perturbations  about  a  solution.  An  obvious  solution  of  this 
set  of  equations  is  that  the  water  surface  is  flat  and  that  the  velocity  is 
zero ,  or 
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H  +  z  =  Constant 


(30) 


U  -  0  (31) 

A  perturbation  about  this  solution  can  be  represented  as 

0  <  u  «  1 

0  <  h  «  1 

where  u  and  h  are  the  perturbation  quantities.  Dropping  all  terms  involv¬ 
ing  products  of  the  perturbation  quantities,  the  resulting  linearized 
equations  are 


uHx  +  Hu*  =  0  (32) 

ghx  -  0  (33) 

The  solution  to  these  equations  is  straightforward.  Certainly  h  is  a  con¬ 
stant,  and  u  can  be  determined  by  supplying  a  depth  variation.  Assume 

H(x)  -  H0£x/L  (34) 

so  that  every  shift  of  a  distance  L  along  x  results  in  the  depth  being 
multiplied  by  £  (the  subscript  o  indicates  the  value  at  the  origin). 
Therefore,  the  solution  is 


u(x)  -  U0£(-x/L>  (35) 

31.  The  question  now  becomes,  what  is  the  numerical  solution  to  this 
same  problem?  Again  using  a  linear  basis  for  Galerkin  finite  elements  and 
assuming  a  numerical  solution  of  the  form: 

uj  -  U0p*  (36) 


the  roots  are  found  to  be 
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p  =  1/C  . 


-(2  +  o 


(37) 


The  first  root  matches  the  analytic  result  at  the  nodes.  The  second  root  is  a 
numerical  artifact;  and  since  it  is  negative,  a  node-to-node  oscillation 
results.  Consider  the  roots  shown  in  Figure  8. 


2  4  .6  .8  1.0  1.2  1.4  1.6  1.8 

Depth  Ratio  per  Element, 


Figure  8.  Depth  sensitivity 

The  upper  curve  is  the  actual  root  and  the  lower  curve  the  spurious  root.  For 
flat  bottom  channels,  £  -  1  and  the  spurious  root  for  £  -  1  is  p  -  -1  . 

If  the  flow  deepens  in  the  direction  of  flow,  corresponding  to  £  >  1  ,  p  is 
less  than  1  and  the  oscillations  will  be  dampened.  However,  if  the  depth 
decreases  In  the  direction  of  flow,  £  <  1  ,  the  oscillations  will  be 
amplified,  possibly  resulting  in  the  model  becoming  unstable.  Therefore, 
flows  in  a  direction  of  decreasing  depth  will  be  prone  to  instability  and  will 
require  additional  resolution  for  the  same  diffusion  coefficient. 
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PART  IV:  CONCLUSIONS 


32.  The  precision  and  stability  of  a  given  numerical  model  depend 
heavily  upon  the  associated  computation  mesh.  When  large  gradients  are  gener¬ 
ated  by  the  inclusion  of  training  structures,  grid-related  problems  are 
magnified.  The  following  guidance  is  noted: 

a.  A  rule  of  thumb  for  the  element  expansion  rate  as  one  moves  away 
from  a  structure  is  L£  =  iLj  .  That  is,  the  element  adjoining 
the  boundary  is  Lj  ,  the  next  element  out  is  of  size  2LX  ,  the 
third  is  3Lj  ,  and  so  on. 

b.  Near  the  dike  tip,  it  is  important  to  orient  the  grid  with  the 
flow  direction  and  make  the  elements  as  nearly  orthogonal  as 
possible . 

c.  Quadrilateral  elements  appear  to  be  more  forgiving  than  trian¬ 
gular  ones  in  cases  in  which  the  orientation  and  skewness  are 
more  extreme. 

d.  Oscillations  in  shallow-water  models  employing  the  mixed  inter¬ 
polation  with  linear  water  surface  and  quadratic  velocity  will 
have  oscillations  if  UL/i/  >  2  .  While  2L  wavelength  oscilla¬ 
tions  in  water  surface  should  not  occur,  oscillations  with 
longer  wavelengths  are  possible,  possibly  leading  to  misleading 
conclusions . 

e.  Depth  decreasing  in  the  direction  of  flow  renders  the  model  less 
stable;  the  converse  is  true  for  depth  increasing. 
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APPENDIX  A:  CONTROL  OF  OSCILLATIONS  IN  TABS-2 


Introduction 


1.  A  common  problem  in  numerical  models  involving  transport  is  the 
presence  of  oscillations  ("wiggles")  in  the  solution.  This  is  a  direct  result 
of  the  numerical  method  as  applied  to  a  first-order  differential  equation  (any 
odd  number  order  actually) .  Often  some  form  of  mixed  interpolation  is  used  to 
alleviate  this  problem  (TABS— 2,  for  example).  In  the  case  of  the  shallow- 
water  equations,  the  water  surface  or  depth  is  interpolated  linearly  and  the 
velocity  is  quadratic.  It  has  been  shown  (Platzman  1978;  Walters  and  Carey 
1983)*  that  this  will  have  a  node-to-node  oscillation  in  the  velocity  field 
but  not  in  the  water  surface;  thus,  including  viscous  effects  in  the  momentum 
equation  can  control  this.  (Note  that  if  the  water  surface  has  oscillations, 
there  is  no  mechanism  to  control  this  with  a  standard  Galerkin  scheme.)  These 
analyses,  however,  did  not  include  viscous  effects  in  the  calculation. 

2.  This  appendix  calculates  the  necessary  viscosity  to  eliminate  oscil¬ 
lations  in  a  mixed  interpolation  finite  element  model  of  the  one-dimensional 
(1-D)  shallow-water  equations  and  the  conditions  under  which  oscillation  in 
water  surface  may  be  present. 


Mixed  Interpolation  Stability  Criteria 


3.  The  discrete  equations  resulting  from  the  Galerkin  weighted-residual 
statement  of  the  1-D  shallow-water  equation  using  mixed  interpolation  are  as 
follows**; 


a.  Conservation  of  mass: 


3U„(-h1  +  h5)  +  h^-Ui  -  4u2  +  4u<,  +  u5)  -  0  (Al) 


*  References  cited  in  this  appendix  are  included  in  the  References  at  the 
end  of  the  main  text. 

**  For  convenience,  symbols  and  unusual  abbreviations  used  in  this  appendix 
are  listed  and  defined  in  the  Notation  (Appendix  B) . 

Al 


b. 


Conservation  of  momentum: 


g(-hi  +  h5)  +  U0(Ui  -  4u2  +  4u4  -  u5) 


+  ^  (2u3  -  16u2  +  28u3  -  16u(,  +  2u5)  »  0  (A2) 


where 


g(— hi  +  h3)  +  ^(-Ui  +  u3)  +  ^(-4^  +  8u2  -  4u3)  =  0  (A3) 

and  a  similar  equation  to  A3  centered  about  node  4 


Uo,hD  =  a  solution  of  the  shallow-water  equations  in  which  UD  and  hc 
are  constants 

h,u  =  perturbation  quantities 
g  =  acceleration  due  to  gravity 
v  =  apparent  viscosity 
d  =  length  of  an  element 

These  equations  are  for  perturbation  equations  about  the  solution 


U  -  U0 

(A4) 

h  =  hQ 

(A5) 

on  the  typical  finite  element  patch,  as  shown  in  Figure  Al. 


|  d  |  d  | 

•  •  ♦  •  • 

i  =  1  2  3  4  5 

Figure  Al .  Finite  element  patch 

All  other  subscripts  indicate  nodal  values.  Here  the  elements  contain  three 
nodes;  u  is  nodally  defined  at  all  three  nodes  while  h  is  defined  only  at 
the  end  nodes  (knots).  Equations  Al  and  A2  result  from  the  weight  function 
for  node  3  applied  against  the  conservation  of  mass  and  momentum  equations, 
respectively.  Equation  A3  is  the  result  of  the  test  function  for  node  2 
applied  to  the  conservation  of  momentum  equation.  A  similar  equation  results 
from  the  application  of  the  test  function  at  node  4;  however,  as  those  ele¬ 
ments  are  uniform,  this  equation  is  identical  to  Equation  3  shifted  one 


A2 


element.  Note  then  that  there  are  three  equations  within  this  patch  enforcing 
momentum  relationships  and  only  one  for  conservation  of  mass.  This  is  a 
direct  consequence  of  the  mixed  interpolation  and  the  use  of  the  Galerkin 
method.  Write  the  nodal  values  as 


=  °p 

(A6) 

(A7) 

A  A 

where  p  is  a  numerical  root,  which  may  be  complex,  and  U  and  h  are 
arbitrary  constants.  Making  these  substitutions  and  finding  the  solutions  to 

A  A 

Equations  A1-A3  for  which  U  =  h  ^  0  yields  the  expressions: 


"  r 

-  . 


|3F2  -  -  1 


9F2  -  21. 


t  ~ 5) *  (*f! '  3tt  ' ']  I5"2*6-?  -5) 


[[' 


r  (3F2  -  6—  ~  1 

f)F2  .  21-E!  -  5 

1-  £f2  -  3.^  +  l]  [3F2  -  6 Z!  -  5] 

LC  R 

r  t  j 

l  T  Jl  R  J 

F2  '  3T  *  L]  ( 

,F*  .  6%  -  1]  - 

3F2  -  «£  -  l]  [>F2  •  3^  •  l]]  - 

(A8) 


P  =  1 


(A9) 


where  F  is  the  Froude  number  and  R  is  the  Peclet  number.  The  second 
expression  (A9)  is  simply  the  statement  that  a  constant  and  h^  is  a 

solution.  Here, 


(A10) 


A3 


(All) 


where  d  is  the  length  of  an  element,  which  corresponds  to  2L  from  para¬ 
graph  27,  main  text.  Equation  A8  then  represents  potential  artificial  solu¬ 
tions.  The  equation  is  resolved  by  realizing  that  most  estuary  and  river 
applications  involve  flow  conditions  characterized  by  small  Froude  numbers. 
Thus  considering  0  <  F  «  1  ,  the  perturbation  equation  is  derived: 


From  this  analysis  one  can  show  that  if  R  =  2  ,  a  node-to— node  oscillation  is 
possible  in  the  velocity  field.  However,  the  water  surface  will  contain  no 
oscillation.  Thus  it  is  desirable  to  have  R  <  2  ,  which  is  essentially  iden¬ 
tical  to  the  well-known  cell  Peclet  criterion  for  linear  basis. 

4.  If  the  cell  Reynolds  number  drops  to  1/2,  an  oscillation  of  wave¬ 
length  2d  is  possible.  Note  that  this  corresponds  to  an  oscillation  in 
water  surface  as  well  as  velocity.  The  amplitude  is  small  for  this  mode,  but 
an  oscillation  in  the  linearized  model  will  exist  in  water  surface. 

Conclusion 
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5.  To  avoid  oscillations  in  TABS-2  model  results,  one  should  make  sure 

that 


R  = 


U.d  „ 
°  <  2 
T7 


(A13) 


This  also  states  that  higher  flow  rates  will  require  increased  eddy  viscosity. 

6.  It  Is  possible  to  have  an  oscillation  in  the  water  surface  even 
without  including  nonlinear  effects. 
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APPENDIX  B :  NOTATION 


C  Concentration 
d  Length  of  an  element 

D  Diffusion  coefficient 

e  Error 
f  Concentration 
fj  Nodal  value  of  f 
F  Froude  number 
g  Acceleration  due  to  gravity 
H  Depth 

j  Node  index  progressing  in  the  x-direction 
L  Distance  between  adjacent  nodes 
R  Peclet  number  or  cell  Reynolds  number 
t  Time 

u,h  Perturbation  quantities 
U  Velocity 

UOI  hQ  Solution  of  the  shallow-water  equations  in  which  UD  and  hQ  are 
constant 

x  Distance  coordinate  direction 

Y  Distance  normal  to  x— axis 

z  Channel  bottom  elevation 
a  Angle 

v  Apparent  viscosity 
p  Numerical  root 

4>  One-dimensional  linear  shape  function 

$  Two-dimensional  linear  basis  functions,  either  triangular  or 
quadrilateral  elements 

Q  Domain 
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